Statistical report of project 2019Apr26_barseqMJ: pairwise comparison(s) of conditions with DESeq2

Table 1: Data files and associated biological conditions.

name	drug	type	plate	new	orig
hom_ctrl_G:A	ctrl	ctrl	G	ctrlG	G1
hom_ctrl_G:B	ctrl	ctrl	G	ctrlG	G2
hom_ctrl_G:C	ctrl	ctrl	G	ctrlG	G3
hom_ctrl_G:D	ctrl	ctrl	G	ctrlG	G4
hom_mtx_G:A	mtx	trt	G	mtx	G5
hom_mtx_G:B	mtx	trt	G	mtx	G6
hom_mtx_G:C	mtx	trt	G	mtx	G7
hom_mtx_G:D	mtx	trt	G	mtx	G8
hom_dox_G:A	dox	trt	G	dox	G9
hom_dox_G:B	dox	trt	G	dox	G10
hom_dox_G:C	dox	trt	G	dox	G11
hom_dox_G:D	dox	trt	G	dox	G12
hom_ctrl_N:A	ctrl	ctrl	N	ctrlN	N1
hom_ctrl_N:B	ctrl	ctrl	N	ctrlN	N2
hom_ctrl_N:C	ctrl	ctrl	N	ctrlN	N3
hom_ma_N:A	ma	trt	N	ma	N4
hom_ma_N:B	ma	trt	N	ma	N5
hom_ma_N:C	ma	trt	N	ma	N6
hom_ndma_N:A	ndma	trt	N	ndma	N7
hom_ndma_N:B	ndma	trt	N	ndma	N8
hom_ndma_N:C	ndma	trt	N	ndma	N9
hom_ndea_N:A	ndea	trt	N	ndea	N10
hom_ndea_N:B	ndea	trt	N	ndea	N11
hom_ndea_N:C	ndea	trt	N	ndea	N12
hom_nmba_N:A	nmba	trt	N	nmba	N13
hom_nmba_N:B	nmba	trt	N	nmba	N14
hom_nmba_N:C	nmba	trt	N	nmba	N15

After loading the data we first have a look at the raw data table itself. The data table contains one row per annotated feature and one column per sequenced sample. Row names of this table are feature IDs (unique identifiers). The table contains raw count values representing the number of reads that map onto the features. For this project, there are 4976 features in the count data table.

Table 2: Partial view of the count data table.

	hom_ctrl_G:A	hom_ctrl_G:B	hom_ctrl_G:C	hom_ctrl_G:D	hom_mtx_G:A	hom_mtx_G:B	hom_mtx_G:C	hom_mtx_G:D	hom_dox_G:A	hom_dox_G:B	hom_dox_G:C	hom_dox_G:D	hom_ctrl_N:A	hom_ctrl_N:B	hom_ctrl_N:C	hom_ma_N:A	hom_ma_N:B	hom_ma_N:C	hom_ndma_N:A	hom_ndma_N:B	hom_ndma_N:C	hom_ndea_N:A	hom_ndea_N:B	hom_ndea_N:C	hom_nmba_N:A	hom_nmba_N:B	hom_nmba_N:C
AAC1	2822	2619	4808	3056	2566	2652	3761	3013	3383	1545	1867	3126	2779	3721	2994	1829	2065	1447	1374	3179	2661	3268	2486	4512	2658	1994	2630
AAC3	0	1	43	0	0	0	0	0	2	1	0	0	0	2	6	0	7	0	8	4	4	9	11	12	5	9	4
AAD14	0	4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
AAD3	295	167	558	246	220	189	300	272	446	147	180	482	267	295	286	132	81	80	168	526	309	443	284	392	111	151	300
AAD4	819	963	1450	541	596	743	1069	678	916	475	417	631	905	938	701	716	651	539	407	1289	863	749	825	1277	752	619	691
AAD6	315	550	704	192	878	472	790	415	833	168	330	346	199	254	305	167	299	184	206	380	309	338	264	580	341	261	91

Looking at the summary of the count table provides a basic description of these raw counts (min and max values, median, etc).

Table 3: Summary of the raw counts.

	hom_ctrl_G:A	hom_ctrl_G:B	hom_ctrl_G:C	hom_ctrl_G:D	hom_mtx_G:A	hom_mtx_G:B	hom_mtx_G:C	hom_mtx_G:D	hom_dox_G:A	hom_dox_G:B	hom_dox_G:C	hom_dox_G:D	hom_ctrl_N:A	hom_ctrl_N:B	hom_ctrl_N:C	hom_ma_N:A	hom_ma_N:B	hom_ma_N:C	hom_ndma_N:A	hom_ndma_N:B	hom_ndma_N:C	hom_ndea_N:A	hom_ndea_N:B	hom_ndea_N:C	hom_nmba_N:A	hom_nmba_N:B	hom_nmba_N:C
Min.	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
1st Qu.	185	200	325	169	236	198	288	178	280	127	162	198	150	211	168	119	104	77	95	235	169	175	125	204	164	123	151
Median	438	452	750	423	495	440	626	407	648	296	365	488	380	530	417	346	305	222	222	571	407	465	362	608	389	301	365
Mean	557	564	943	545	604	553	773	524	785	357	441	606	501	691	543	500	447	326	284	747	537	606	488	827	493	388	478
3rd Qu.	793	788	1317	758	832	766	1065	727	1101	496	609	843	704	963	743	709	629	461	390	1032	749	844	686	1162	682	540	661
Max.	7327	7754	13266	7327	8876	7499	10583	7640	30631	13985	17098	21588	7061	9694	8320	7864	6531	5758	3751	10546	7570	11535	9546	15697	6288	5747	6924

Figure 1 shows the total number of mapped and counted reads for each sample. We expect total read counts to be similar within conditions, they may be different across conditions. Total counts sometimes vary widely between replicates. This may happen for several reasons, including:

• different rRNA contamination levels between samples (even between biological replicates);
• slight differences between library concentrations, since they may be difficult to measure with high precision.

Figure 1: Number of mapped reads per sample. Colors refer to the biological condition of the sample.

Figure 2 shows the proportion of features with no read count in each sample. We expect this proportion to be similar within conditions. Features with null read counts in the 27 samples are left in the data but are not taken into account for the analysis with DESeq2. Here, 0 features (0%) are in this situation (dashed line). Results for those features (fold-change and p-values) are set to NA in the results files.

Figure 2: Proportion of features with null read counts in each sample.

Figure 3 shows the distribution of read counts for each sample. For sake of readability, \(\text{log}_2(\text{counts}+1)\) are used instead of raw counts. Again we expect replicates to have similar distributions. In addition, this figure shows if read counts are preferably low, medium or high. This depends on the organisms as well as the biological conditions under consideration.

Figure 3a: Density distribution of read counts.

Figure 3b: Density distribution of raw and normalized counts.

It may happen that one or a few features capture a high proportion of reads (up to 20% or more). This phenomenon should not influence the normalization process. The DESeq2 normalization has proved to be robust to this situation [Dillies, 2012]. Anyway, we expect these high count features to be the same across replicates. They are not necessarily the same across conditions. Figure 4 and table 4 illustrate the possible presence of such high count features in the data set. - note: rule of thumb, anything above 0.2 is a likely outlier

Figure 4: Percentage of reads associated with the sequence having the highest count (provided in each box on the graph) for each sample.

Table 4: Percentage of reads associated with the sequences having the highest counts.

	hom_ctrl_G:A	hom_ctrl_G:B	hom_ctrl_G:C	hom_ctrl_G:D	hom_mtx_G:A	hom_mtx_G:B	hom_mtx_G:C	hom_mtx_G:D	hom_dox_G:A	hom_dox_G:B	hom_dox_G:C	hom_dox_G:D	hom_ctrl_N:A	hom_ctrl_N:B	hom_ctrl_N:C	hom_ma_N:A	hom_ma_N:B	hom_ma_N:C	hom_ndma_N:A	hom_ndma_N:B	hom_ndma_N:C	hom_ndea_N:A	hom_ndea_N:B	hom_ndea_N:C	hom_nmba_N:A	hom_nmba_N:B	hom_nmba_N:C
COX5A	0.26	0.28	0.28	0.27	0.26	0.27	0.26	0.29	0.22	0.24	0.23	0.26	0.28	0.28	0.31	0.32	0.29	0.35	0.27	0.27	0.28	0.35	0.36	0.33	0.26	0.30	0.29
YDR417C	0.21	0.22	0.20	0.20	0.22	0.22	0.23	0.21	0.09	0.09	0.11	0.12	0.20	0.21	0.20	0.16	0.16	0.18	0.16	0.15	0.19	0.12	0.11	0.12	0.18	0.18	0.19
KGD1	0.20	0.25	0.23	0.23	0.29	0.23	0.24	0.21	0.22	0.27	0.22	0.23	0.21	0.23	0.23	0.25	0.24	0.24	0.25	0.28	0.26	0.22	0.16	0.17	0.22	0.23	0.22
EMI1|48	0.19	0.17	0.18	0.21	0.29	0.26	0.28	0.24	0.16	0.19	0.17	0.21	0.25	0.25	0.26	0.20	0.22	0.23	0.26	0.23	0.21	0.28	0.29	0.29	0.24	0.23	0.26
RTC6|0018	0.16	0.16	0.16	0.17	0.25	0.21	0.25	0.20	0.23	0.23	0.24	0.22	0.18	0.16	0.18	0.13	0.15	0.17	0.19	0.20	0.18	0.12	0.11	0.09	0.17	0.16	0.18
SSN2	0.17	0.16	0.18	0.16	0.16	0.18	0.17	0.16	0.78	0.79	0.78	0.71	0.15	0.17	0.17	0.31	0.28	0.24	0.21	0.18	0.19	0.38	0.39	0.38	0.18	0.17	0.18
YDR442W	0.06	0.07	0.07	0.06	0.06	0.05	0.06	0.06	0.23	0.20	0.18	0.19	0.06	0.05	0.05	0.07	0.07	0.07	0.05	0.05	0.04	0.09	0.07	0.07	0.07	0.05	0.05

We may wish to assess the similarity between samples across conditions. A pairwise scatter plot is produced (figure 5) to show how replicates and samples from different biological conditions are similar or different (\(\text{log}_2(\text{counts}+1)\) are used instead of raw count values). Moreover, as the Pearson correlation has been shown not to be relevant to measure the similarity between replicates, the SERE statistic has been proposed as a similarity index between RNA-Seq samples [1]. It measures whether the variability between samples is random Poisson variability or higher. Pairwise SERE values are printed in the lower triangle of the pairwise scatter plot. The value of the SERE statistic is:

• 0 when samples are identical (no variability at all: this may happen in the case of a sample duplication);

• 1 for technical replicates (technical variability follows a Poisson distribution);

• greater than 1 for biological replicates and samples from different biological conditions (biological variability is higher than technical one, data are over-dispersed with respect to Poisson). The higher the SERE value, the lower the similarity. It is expected to be lower between biological replicates than between samples of different biological conditions. Hence, the SERE statistic can be used to detect inversions between samples.

Figure 5: Pairwise comparison of samples (not produced when more than 30 samples).

Variability within the experiment: data exploration

The main variability within the experiment is expected to come from biological differences between the samples. This can be checked in two ways. The first one is to perform a hierarchical clustering of the whole sample set. This is performed after a transformation of the count data which can be either a Variance Stabilizing Transformation (VST) or a regularized log transformation (rlog) [2,3].

A VST is a transformation of the data that makes them homoscedastic, meaning that the variance is then independent of the mean. It is performed in two steps: (i) a mean-variance relationship is estimated from the data with the same function that is used to normalize count data and (ii) from this relationship, a transformation of the data is performed in order to get a dataset in which the variance is independent of the mean. The homoscedasticity is a prerequisite for the use of some data analysis methods, such as hierarchical clustering or Principal Component Analysis (PCA). The regularized log transformation is based on a GLM (Generalized Linear Model) on the counts and has the same goal as a VST but is more robust in the case when the size factors vary widely.

Figure 6 shows the dendrogram obtained from log transformed data. An euclidean distance is computed between samples, and the dendrogram is built upon the Ward criterion. We expect this dendrogram to group replicates and separate biological conditions.

Figure 6: Sample clustering based on normalized data.

Another way of visualizing the experiment variability is to look at the first principal components of the PCA, as shown on the figure 7. On this figure, the first principal component (PC1) is expected to separate samples from the different biological conditions, meaning that the biological variability is the main source of variance in the data.

Figure 7: First two components of a Principal Component Analysis, with percentages of variance associated with each axis.

Another way of visualizing the experiment variability is to look at the first two dimensions of a multidimensional scaling plot, as shown on figure 7. On this figure, the first dimension is expected to separate samples from the different biological conditions, meaning that the biological variability is the main source of variance in the data.

Figure 7b: Multidimensional scaling plot of the samples.

Differential analysis

Modelisation

DESeq2 aims at fitting one linear model per feature. For this project, the design used is counts ~ drug and the goal is to estimate the models’ coefficients which can be interpreted as \(\log_2(\texttt{FC})\). These coefficients will then be tested to get p-values and adjusted p-values.

Outlier detection

Model outliers are features for which at least one sample seems unrelated to the experimental or study design. For every feature and for every sample, the Cook’s distance [4] reflects how the sample matches the model. A large value of the Cook’s distance indicates an outlier count and p-values are not computed for the corresponding feature.

Dispersions estimation

The DESeq2 model assumes that the count data follow a negative binomial distribution which is a robust alternative to the Poisson law when data are over-dispersed (the variance is higher than the mean). The first step of the statistical procedure is to estimate the dispersion of the data. Its purpose is to determine the shape of the mean-variance relationship. The default is to apply a GLM (Generalized Linear Model) based method (fitType=“parametric”), which can handle complex designs but may not converge in some cases. The alternative is to use fitType=“local” as described in the original paper [2]. The parameter used for this project is fitType=“local”. Then, DESeq2 imposes a Cox Reid-adjusted profile likelihood maximization [5] and uses the maximum a posteriori (MAP) of the dispersion [Wu, 2013].

Figure 8: Diagnostic of the estimation of the size factors.

The figure 9 shows that the scaling factors of DESeq2 and the total count normalization factors may not perform similarly.

Figure 9: Plot of the estimated size factors and the total number of reads per sample.

Boxplots are often used as a qualitative measure of the quality of the normalization process, as they show how distributions are globally affected during this process. We expect normalization to stabilize distributions across samples. Figure 10 shows boxplots of raw (left) and normalized (right) data respectively.

Figure 10: Boxplots of raw (left) and normalized (right) read counts.

Differential analysis

Modelisation

DESeq2 aims at fitting one linear model per feature. For this project, the design used is counts ~ drug and the goal is to estimate the models’ coefficients which can be interpreted as \(\log_2(\texttt{FC})\). These coefficients will then be tested to get p-values and adjusted p-values.

Outlier detection

Model outliers are features for which at least one sample seems unrelated to the experimental or study design. For every feature and for every sample, the Cook’s distance [4] reflects how the sample matches the model. A large value of the Cook’s distance indicates an outlier count and p-values are not computed for the corresponding feature.

Dispersions estimation

The DESeq2 model assumes that the count data follow a negative binomial distribution which is a robust alternative to the Poisson law when data are over-dispersed (the variance is higher than the mean). The first step of the statistical procedure is to estimate the dispersion of the data. Its purpose is to determine the shape of the mean-variance relationship. The default is to apply a GLM (Generalized Linear Model) based method (fitType=“parametric”), which can handle complex designs but may not converge in some cases. The alternative is to use fitType=“local” as described in the original paper [2]. The parameter used for this project is fitType= local. Then, DESeq2 imposes a Cox Reid-adjusted profile likelihood maximization [5] and uses the maximum (MAP) of the dispersion [Wu, 2013].

Figure 11: Dispersion estimates (left) and diagnostic of log-normality (right).

The left panel on figure 11 shows the result of the dispersion estimation step. The x- and y-axes represent the mean count value and the estimated dispersion respectively. Black dots represent empirical dispersion estimates for each feature (from the observed counts). The red dots show the mean-variance relationship function (fitted dispersion value) as estimated by the model. The blue dots are the final estimates from the maximum a posteriori and are used to perform the statistical test. Blue circles (if any) point out dispersion outliers. These are features with a very high empirical variance (computed from observed counts). These high dispersion values fall far from the model estimation. For these features, the statistical test is based on the empirical variance in order to be more conservative than with the MAP dispersion. These features will have low chance to be declared significant. The figure on the right panel allows to check the hypothesis of log-normality of the dispersions.

Statistical test for differential expression

Once the dispersion estimation and the model fitting have been done, DESeq2 can perform the statistical testing. Figure 12 shows the distributions of raw p-values computed by the statistical test for the comparison(s) done. This distribution is expected to be a mixture of a uniform distribution on \([0,1]\) and a peak around 0 corresponding to the differentially expressed features.

Figure 12: Distribution(s) of raw p-values.

Independent filtering

DESeq2 can perform an independent filtering to increase the detection power of differentially expressed features at the same experiment-wide type I error. Since features with very low counts are not likely to see significant differences typically due to high dispersion, it defines a threshold on the mean of the normalized counts irrespective of the biological condition. This procedure is independent because the information about the variables in the design formula is not used [3].

Final results

A p-value adjustment is performed to take into account multiple testing and control the false positive rate to a chosen level \(\alpha\). For this analysis, a BH p-value adjustment was performed [6] and the level of controlled false positive rate was set to 0.05.

Figure 13 are the fitness defects (\(\log_2\) ratios) plotted as a function of mean counts – this is useful for knowing whether significant are coming from unusual low (noisy) counts

Figure 13: Log Ratio vs. mean CPM of each comparison.

Figure 14 shows the volcano plots for the comparisons performed and differentially expressed features are still highlighted in red. A volcano plot represents the log of the adjusted P value as a function of the log ratio of differential expression.

Figure 14: Volcano plot(s) of each comparison. Red dots represent significantly differentially expressed features.

Full results as well as lists of differentially expressed features are provided in the following text files which can be easily read in a spreadsheet. For each comparison:

• TestVsRef.complete.txt contains results for all the features;
• TestVsRef.up.txt contains results for significantly up-regulated features. Features are ordered from the most significant adjusted p-value to the less significant one;
• TestVsRef.down.txt contains results for significantly down-regulated features. Features are ordered from the most significant adjusted p-value to the less significant one.

These files contain the following columns:

• Id: unique feature identifier;
• sampleName: raw counts per sample;
• norm.sampleName: rounded normalized counts per sample;
• baseMean: base mean over all samples;
• ctrl, dox, ma, mtx, ndea, ndma and nmba: means (rounded) of normalized counts of the biological conditions;
• FoldChange: fold change of expression, calculated as \(2^{\log_2(\text{FC})}\);
• log2FoldChange: \(\log_2(\text{FC})\) as estimated by the GLM model. It reflects the differential expression between Test and Ref and can be interpreted as \(\log_2(\frac{\text{Test}}{\text{Ref}})\). If this value is:
  • around 0: the feature expression is similar in both conditions;
  • positive: the feature is up-regulated (\(\text{Test} > \text{Ref}\));
  • negative: the feature is down-regulated (\(\text{Test} < \text{Ref}\));
• pvalue: raw p-value from the statistical test;
• padj: adjusted p-value on which the cut-off \(\alpha\) is applied;
• dispGeneEst: dispersion parameter estimated from feature counts (i.e. black dots on figure 11);
• dispFit: dispersion parameter estimated from the model (i.e. red dots on figure 11);
• dispMAP: dispersion parameter estimated from the Maximum A Posteriori model;
• dispersion: final dispersion parameter used to perform the test (i.e. blue dots and circles on figure 11);
• betaConv: convergence of the coefficients of the model (TRUE or FALSE);
• maxCooks: maximum Cook’s distance of the feature.

R session information and parameters

The versions of the R software and Bioconductor packages used for this analysis are listed below. It is important to save them if one wants to re-perform the analysis in the same conditions.

• R version 3.5.2 (2018-12-20), x86_64-redhat-linux-gnu
  
• Locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8, LC_IDENTIFICATION=C
  
• Running under: CentOS release 6.10 (Final)
  
• Matrix products: default
  
• BLAS: /usr/lib64/R/lib/libRblas.so
  
• LAPACK: /usr/lib64/R/lib/libRlapack.so
  
• Base packages: base, datasets, graphics, grDevices, grid, methods, parallel, stats, stats4, utils
  
• Other packages: Biobase 2.42.0, BiocGenerics 0.28.0, BiocParallel 1.16.6, DelayedArray 0.8.0, dendextend 1.9.0, DESeq2 1.22.2, dplyr 0.8.0.1, edgeR 3.24.3, GenomeInfoDb 1.18.2, GenomicRanges 1.34.0, ggplot2 3.1.0, ggrepel 0.8.0, IRanges 2.16.0, limma 3.38.3, matrixStats 0.54.0, S4Vectors 0.20.1, SARTools 1.6.6, SummarizedExperiment 1.12.0, xtable 1.8-4
  
• Loaded via a namespace (and not attached): acepack 1.4.1, annotate 1.60.1, AnnotationDbi 1.44.0, assertthat 0.2.1, backports 1.1.4, base64enc 0.1-3, bit 1.1-14, bit64 0.9-7, bitops 1.0-6, blob 1.1.1, checkmate 1.9.1, class 7.3-15, cluster 2.0.8, colorspace 1.4-1, compiler 3.5.2, crayon 1.3.4, data.table 1.12.2, DBI 1.0.0, DEoptimR 1.0-8, digest 0.6.18, diptest 0.75-7, evaluate 0.13, flexmix 2.3-15, foreign 0.8-71, Formula 1.2-3, fpc 2.1-11.1, genefilter 1.64.0, geneplotter 1.60.0, GenomeInfoDbData 1.2.0, glue 1.3.1, gridExtra 2.3, gtable 0.3.0, highr 0.8, Hmisc 4.2-0, htmlTable 1.13.1, htmltools 0.3.6, htmlwidgets 1.3, kernlab 0.9-27, knitr 1.22, labeling 0.3, lattice 0.20-38, latticeExtra 0.6-28, lazyeval 0.2.2, locfit 1.5-9.1, magrittr 1.5, MASS 7.3-51.4, Matrix 1.2-17, mclust 5.4.2, memoise 1.1.0, modeltools 0.2-22, munsell 0.5.0, mvtnorm 1.0-10, nnet 7.3-12, pillar 1.3.1, pkgconfig 2.0.2, plyr 1.8.4, prabclus 2.2-7, purrr 0.3.1, R6 2.4.0, RColorBrewer 1.1-2, Rcpp 1.0.1, RCurl 1.95-4.12, rlang 0.3.1, rmarkdown 1.12, robustbase 0.93-3, rpart 4.1-15, RSQLite 2.1.1, rstudioapi 0.9.0, scales 1.0.0, splines 3.5.2, stringi 1.4.3, stringr 1.4.0, survival 2.44-1.1, tibble 2.0.1, tidyselect 0.2.5, tools 3.5.2, trimcluster 0.1-2.1, viridis 0.5.1, viridisLite 0.3.0, whisker 0.3-2, withr 2.1.2, xfun 0.6, XML 3.98-1.19, XVector 0.22.0, yaml 2.2.0, zlibbioc 1.28.0

Parameter values used for this analysis are:

• workDir: /home/ggiaever/RProjects/2019barseq_marjan
• projectName: 2019Apr26_barseqMJ
• author: gg
• targetFile: /home/common/barseq_gg_pdata_exp_info/2019april20_mar29_dox_ndea_pdata.txt
• rawDir: /home/common/barseq_output
• featuresToRemove: NULL
• varInt: drug
• condRef: ctrl
• batch: NULL
• fitType: local
• cooksCutoff: TRUE
• independentFiltering: TRUE
• alpha: 0.05
• pAdjustMethod: BH
• typeTrans: VST
• locfunc: median
• colors: darkorange, dodgerblue, limegreen, navy, mediumpurple, royalblue3, darkolivegreen4, firebrick, cyan, hotpink, plum4, blue, magenta2, skyblue, green, red, steelblue, tomato, purple, yellow

Bibliography

1. Schulze SK, Kanwar R, Gölzenleuchter M et al. SERE: Single-parameter quality control and sample comparison for rna-seq. BMC Genomics 2012 ; 13: 524.

2. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biology 2010 ; 11: R106.

3. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for rna-seq data with deseq2. Genome Biology 2014 ; 15: 550.

4. Cook RD. Detection of influential observation in linear regression. Technometrics 1977 ; 19: 15–18.

5. Cox DR, Reid N. Parameter orthogonality and approximate conditional inference. Journal of the Royal Statistical Society. Series B (Methodological) 1987 ; 49: 1–39.

6. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. Series B (Methodological) 1995 ; 57: 289–300.